1 C PSIS PROGRAM
2 DIMENSION TN(6,10,10),AP(10),PN(7),FEXP(10),PSIR(10),PSITH(10),
3 1PSIFH(10),X(7),Y(7)
4 DOUBLE PRECISION BHETA,PN,DIST,THETA,FHI,X,Y,FCT,DIFSC,RAD,AMPL,
5 1BHETAD,FHID,THETAD,PI,KV,AMPL2,A,B,C,D,E,deltheta
6 COMPLEX*16 TN,AP,PSIR,PSITH,PSIFH,FEXP,Q1,Q2,P1,P2,R1,R2,FC
7 COMMON/FAC/FCT(100),NRISK,NTRIX
8 NAMELIST/HEJ/Q1,Q2,P1,P2,R1,R2
9 c CALL UNSTAE
10 PI = 4.0D0*DATAN(1.0D0)
11 NRISK = 57
12 NTRIX = 39
13 NEND = 78
14
15 c polarization index 0=orthogonal 1=parallel
16 NP = 1
17
18 c parameter = 1 if other paralization is needed
19 c NPCHAN = 1
20 NPCHAN = 1
21
22 c dimension of t-matrix
23 ND = 10
24
25 c parameter for VPSI for computation of basis functions
26 NB = 1
27
28 c wave vector
29 KV = 2.0D0*PI/15.53333D0
30
31 c polar angle (rad) of incident wave
32 BHETA = PI/2.0D0
33 BHETAD = BHETA*180.0D0/PI
34
35 c wave vector time radius to observation point
36 DIST = 1.0D0
37
38
39 c azimut angle (rad)
40 FHI = 0.0D0
41 c azimut angle (degrees)
42 FHID = FHI*180.0D0/PI
43
44 N = NRISK
45 L = NTRIX
46 M = NEND-2
47 N1 = N-1
48 DO 5 K = 1,100
49 5 FCT(K) = 0.0D0
50 FCT(1) = 1.0D0
51 DO 6 K = 1,N1
52 6 FCT(K+1) = FCT(K)*DFLOAT(K)
53 FCT(N+1) = 0.1D0**L*FCT(N)*DFLOAT(N)
54 DO 7 K = N,M
55 7 FCT (K+2) = FCT(K+1)*DFLOAT(K+1)
56 NPC = NPCHAN+1
57 GO TO (11,12), NPC
58 11 CONTINUE
59 FC = (1.0D0,0.0D0)
60 GO TO 13
61 12 CONTINUE
62 FC = (-1.0D0,0.0D0)
63 13 CONTINUE
64 READ (11) TN
65
66 ccc
67 deltheta=PI/180.0
68 theta=0.0
69 do 100 ith=1,180
70 ith1=ith-1
71 theta = theta + deltheta
72
73 c polar angle (rad)
74 c THETA = PI/2.0D0
75 c polar angle (degrees)
76 THETAD = THETA*180.0D0/PI
77 ccc
78
79
80 NS = ND/2+1
81 P1 = (0.0D0,0.0D0)
82 P2 = (0.0D0,0.0D0)
83 DO 23 M1 = 1,NS
84 M = M1-1
85 MD = 2*M-1
86 IF((M-2).LT.0) MD = 1
87 CALL VKOEF(BHETA,NP,M,ND,AP,PN)
88 CALL VPSI(DIST,THETA,FHI,NP,NB,M,ND,PSIR,PSITH,PSIFH,PN,X,Y)
89 DO 21 I = MD,ND
90 Q1 = (0.0D0,0.0D0)
91 DO 20 J = MD,ND
92 Q1 = Q1+TN(M1,I,J)*AP(J)*FC**(I+J)
93 C Seite 20
94
95
96 20 CONTINUE
97
98
99 FEXP(I) = Q1
100 21 CONTINUE
101 Q1 = (0.0D0,0.0D0)
102 Q2 = (0.0D0,0.0D0)
103 DO 22 I = MD,ND
104 Q1 = Q1+PSITH(I)*FEXP(I)
105 Q2 = Q2+PSIFH(I)*FEXP(I)
106 22 CONTINUE
107 P1 = P1+Q1
108 P2 = P2+Q2
109 R1 = P1/DCMPLX(KV,0.0D0)
110 R2 = P2/DCMPLX(KV,0.0D0)
111 WRITE(6,HEJ)
112 AMPL2= CDABS(P1)**2+CDABS(P2)**2
113 AMPL = DSQRT(AMPL2)
114 A = THETAD
115 B = FHID
116 C = AMPL
117 D = AMPL2
118 E = AMPL/KV
119 c PRINT 40,A,B,C,D,E
120 c 40 FORMAT(5D25.15)
121 PRINT 40,A
122 40 FORMAT(5D17.7)
123
124 23 CONTINUE
125
126 write (12, 50) a, d
127
128 50 format(D17.7,3X,D17.7)
129 100 CONTINUE
130 300 CONTINUE
131 999 CONTINUE
132 STOP
133 cc DEBUG SUBCHK
134 END
135 C Seite 21
136
137
138
139 SUBROUTINE BESSEL(NORDER,ARGMNT,ANSWR,IERROR)
140 DOUBLE PRECISION ARGMNT,ANSWR,X,SUM,APR,TOPR,CI,CNI,ACR,PROD,FACT
141 IERROR = 0
142 N = NORDER
143 X = ARGMNT
144 SUM = 1.0D0
145 APR = 1.0D0
146 TOPR = -0.5D0*X*X
147 CI = 1.0D0
148 CNI = DFLOAT(2*N+3)
149 DO 60 I = 1,100
150 ACR = TOPR*APR/(CI*CNI)
151 SUM = SUM+ACR
152 IF(DABS(ACR/SUM)-1.0D-20)100,100,40
153 40 APR = ACR
154 CI = CI+1.0D0
155 CNI = CNI+2.0D0
156 60 CONTINUE
157 IERROR = I
158 PRINT 10
159 10 FORMAT(24HO ERROR IN SUM OF BESSEL)
160 GO TO 200
161 100 PROD = DFLOAT(2*N+1)
162 FACT = 1.0D0
163 IF(N)160,160,120
164 120 DO 140 IFCT = 1,N
165 FACT = FACT*X/PROD
166 PROD = PROD-2.0D0
167 140 CONTINUE
168 160 ANSWR = FACT*SUM
169 200 RETURN
170 END
171 C Seite 23
172
173
174 SUBROUTINE BN(PCKR,NRINK,BSSLSP,CNEUMN)
175 DIMENSION BSSLSP(NRINK),CNEUMN(NRINK)
176 DOUBLE PRECISION PCKR,ANSWR,ANSA,ANSB,ANSC,CONN,CMULN,SNSA,SNSB,
177 1SNSC,BSSLSP,CNEUMN
178 NRANKI = NRINK
179 NRANK = NRANKI-1
180 NVAL = NRANK-1
181 DO 40 I = 1,4
182 CALL BESSEL(NVAL,PCKR,ANSWR,IERROR)
183 IF(IERROR)20,20,32
184 20 ANSA = ANSWR
185 NVAL = NVAL+1
186 CALL BESSEL(NVAL,PCKR,ANSWR,IERROR)
187 IF(IERROR)24,24,28
188 24 ANSB = ANSWR
189 GO TO 60
190 28 NVAL = NVAL-1
191 32 NVAL = NVAL+NRANK
192 40 CONTINUE
193 STOP
194 60 IF(NVAL-NRANK)100,100,64
195 64 IEND = NVAL-NRANK
196 CONN = DFLOAT(2*(NVAL-1)+1)
197 DO 72 IP = 1,IEND
198 ANSC = CONN*ANSA/PCKR-ANSB
199 CONN = CONN-2.0D0
200 ANSB = ANSA
201 ANSA = ANSC
202 72 CONTINUE
203 100 BSSLSP(NRANKI) = ANSB
204 BSSLSP(NRANKI-1) = ANSA
205 CONN = DFLOAT(NRANK+NRANK-1)
206 IE = NRANKI-2
207 JE = IE
208 DO 120 JB = 1,JE
209 ANSC = CONN*ANSA/PCKR-ANSB
210 BSSLSP(IE) = ANSC
211 ANSB = ANSA
212 ANSA = ANSC
213 IE = IE-1
214 CONN = CONN-2.0D0
215 120 CONTINUE
216 CMULN = 3.0D0
217 SNSA = -DCOS(PCKR)/PCKR
218 SNSB = -DCOS(PCKR)/(PCKR*PCKR)-DSIN(PCKR)/PCKR
219 CNEUMN(1) = SNSA
220 CNEUMN(2) = SNSB
221 DO 280 I = 3,NRANKI
222 SNSC = CMULN*SNSB/PCKR-SNSA
223 CNEUMN(I) = SNSC
224 SNSA = SNSB
225 SNSB = SNSC
226 CMULN = CMULN+2.0D0
227 280 CONTINUE
228 RETURN
229 END
230 C Seite 24
231
232
233 SUBROUTINE CBESS (NORDER,ARGMNT,ANSWR,IERROR)
234 COMPLEX*16 ARGMNT,ANSWR,X,SUM,APR,TOPR,CI,CNI,ACR,PROD,FACT
235 IERROR = 0
236 N = NORDER
237 X = ARGMNT
238 SUM = (1.0D0,0.0D0)
239 APR = (1.0D0,0.0D0)
240 TOPR = -(0.5D0,0.0D0)*X*X
241 CI = (1.0D0,0.0D0)
242 CNI = DCMPLX(DFLOAT(2*N+3),0.0D0)
243 DO 60 I = 1,100
244 ACR = TOPR*APR/(CI*CNI)
245 SUM = SUM+ACR
246 IF(CDABS(ACR/SUM)-1.0D-20)100,100,40
247 40 APR = ACR
248 CI = CI+(1.0D0,0.0D0)
249 CNI = CNI+(2.0D0,0.0D0)
250 60 CONTINUE
251 IERROR = 1
252 PRINT 10
253 10 FORMAT('0 ERROR IN SUM OF CBESS')
254 GO TO 200
255 100 PROD = DCMPLX(DFLOAT(2*N+1),0.0D0)
256 FACT = (1.0D0,0.0D0)
257 IF(N)160,160,120
258 120 DO 140 IFCT = 1,N
259 FACT = FACT*X/PROD
260 PROD = PROD-(2.0D0,0.0D0)
261 140 CONTINUE
262 160 ANSWR = FACT*SUM
263 200 RETURN
264 END
265 C Seite 25
266
267
268 SUBROUTINE CBN(PCKR,NRINK,BSSLSP,CNEUMN)
269 DIMENSION BSSLSP(NRINK),CNEUMN(NRINK)
270 COMPLEX*16 PCKR,ANSWR,ANSA,ANSB,ANSC,CONN,CMULN,SNSA,SNSB,SNSC,
271 1BSSLSP,CNEUMN
272 NRANKI = NRINK
273 NRANK = NRANKI-1
274 NVAL = NRANK-1
275 DO 40 I = 1,4
276 CALL CBESS(NVAL,PCKR,ANSWR,IERROR)
277 IF(IERROR)20,20,32
278 20 ANSA = ANSWR
279 NVAL = NVAL+1
280 CALL CBESS(NVAL,PCKR,ANSWR,IERROR)
281 IF(IERROR)24,24,28
282 24 ANSB = ANSWR
283 GO TO 60
284 28 NVAL = NVAL-1
285 32 NVAL = NVAL+NRANK
286 40 CONTINUE
287 STOP
288 60 IF(NVAL-NRANK)100,100,64
289 64 IEND = NVAL-NRANK
290 CONN = DCMPLX(DFLOAT(2*(NVAL-1)+1),0.0D0)
291 DO 72 IP = 1,IEND
292 ANSC = CONN*ANSA/PCKR-ANSB
293 CONN = CONN-(2.0D0,0.0D0)
294 ANSB = ANSA
295 ANSA = ANSC
296 72 CONTINUE
297 100 BSSLSP(NRANKI) = ANSB
298 BSSLSP(NRANKI-1) = ANSA
299 CONN = DCMPLX(DFLOAT(NRANK+NRANK-1),0.0D0)
300 IE = NRANKI-2
301 JE = IE
302 DO 120 JB = 1,JE
303 ANSC = CONN*ANSA/PCKR-ANSB
304 BSSLSP(IE) = ANSC
305 ANSB = ANSA
306 ANSA = ANSC
307 IE = IE-1
308 CONN = CONN-(2.0D0,0.0D0)
309 120 CONTINUE
310 CMULN = (3.0D0,0.0D0)
311 SNSA = -CDCOS(PCKR)/PCKR
312 SNSB = -CDCOS(PCKR)/(PCKR*PCKR)-CDSIN(PCKR)/PCKR
313 CNEUMN(1) = SNSA
314 CNEUMN(2) = SNSB
315 DO 280 I = 3,NRANKI
316 SNSC = CMULN*SNSB/PCKR-SNSA
317 CNEUMN(I) = SNSC
318 SNSA = SNSB
319 SNSB = SNSC
320 CMULN = CMULN+(2.0D0,0.0D0)
321 280 CONTINUE
322 RETURN
323 END
324 C Seite 26
325
326
327 SUBROUTINE LEG(THETA,M,NRJNK,PNMLLG)
328 DIMENSION PNMLLG(NRJNK)
329 DOUBLE PRECISION THETA,PLA,PLB,PLC,DTWM,CNM,CNMM,CNMUL,PRODM,
330 1QUANM,PNMLLG
331 NRANKI = NRJNK
332 KMV = M
333 DTWM = DFLOAT (2*M+1)
334 IF(THETA)16,4,16
335 4 IF(KMV)12,12,6
336 6 DO 8 ILG = 1,NRANKI
337 PNMLLG(ILG) = 0.0D0
338 8 CONTINUE
339 GO TO 88
340 12 PNMLLG(1) = 1.0D0
341 PLA = 1.0D0
342 GO TO 48
343 16 IF(KMV)20,20,40
344 20 PLA = 1.0D0
345 PLB = DCOS(THETA)*PLA
346 PNMLLG(1) = PLA
347 PNMLLG(2) = PLB
348 IBEG = 3
349 GO TO 60
350 40 DO 44 ILG = 1,KMV
351 PNMLLG(ILG) = 0.0D0
352 44 CONTINUE
353 PRODM = 1.0D0
354 QUANM = DFLOAT(KMV)
355 DO 52 IFCT = 1,KMV
356 QUANM = QUANM+1.0D0
357 PRODM = QUANM*PRODM/2.0D0
358 52 CONTINUE
359 PLA = PRODM*DSIN(THETA)**KMV
360 PNMLLG(KMV+1) = PLA
361 48 PLB = DTWM*DCOS(THETA)*PLA
362 PNMLLG(KMV+2) = PLB
363 IBEG = KMV+3
364 60 CNMUL = DFLOAT(IBEG+IBEG-3)
365 CNM = 2.0D0
366 CNMM = DTWM
367 DO 80 ILGR = IBEG,NRANKI
368 PLC = (CNMUL*DCOS(THETA)*PLB-CNMM*PLA)/CNM
369 PNMLLG(ILGR) = PLC
370 PLA = PLB
371 PLB = PLC
372 CNMUL = CNMUL+2.0D0
373 CNM = CNM+1.0D0
374 CNMM = CNMM+1.0D0
375 80 CONTINUE
376 88 RETURN
377 END
378 C Seite 27
379
380
381 FUNCTION TRIXJ(J1,J2,J,M,FCT,N,L)
382 IMPLICIT REAL*8(A-H,O-Z)
383 DIMENSION FCT(1)
384 INTEGER Z,ZMIN,ZMAX
385 Y=1.0D0
386 CC=0.0D0
387 JSUM=J1+J2+J
388 JM1=J1-IABS(M)
389 JM2=J2-IABS(M)
390 IF((MOD(JSUM,2).NE.0).OR.(MOD(JM1,2).NE.0).OR.(MOD(JM2,2).NE.0)
391 1.OR.(JM1.LT.0).OR.(JM2.LT.0))
392 2GO TO 3
393 IF((J.GT.J1+J2).OR.(L.LT.IABS(J1-J2))) GO TO 1
394 ZMIN=0
395 IF(J-J2+M.LT.0) ZMIN=-J+J2-M
396 IF(J-J1+M+ZMIN.LT.0) ZMIN=-J+J1-M
397 ZMAX=J1+J2-J
398 IF(J2-M-ZMAX.LT.0) ZMAX=J2-M
399 IF(J1-M-ZMAX.LT.0) ZMAX=J1-M
400 JA=(J1+M)/2+1
401 JB=JA-M
402 JC=(J2-M)/2+1
403 JD=JC+M
404 JE=J/2+1
405 JF=(J1+J2-J)/2+1
406 JG=JA+JB-JF
407 JH=JC+JD-JF
408 JJ=2*JE+JF-1
409 IF(JJ.GT.N) Y=0.1D0**L
410 IF(FCT(JJ))7,5,7
411 7 CONTINUE
412 IA=ZMIN/2
413 IB=JF-IA+1
414 IC=JB-IA+1
415 ID=JC-IA+1
416 IE=JA-JF+IA
417 IF=JD-JF+IA
418 FASE=1.0D0
419 IF(MOD(IA,2).EQ.0) FASE=-FASE
420 Z=ZMIN
421 2 IA=IA+1
422 IB=IB-1
423 IC=IC-1
424 ID=ID-1
425 IE=IE+1
426 IF=IF+1
427 FASE=-FASE
428 CC=CC+FASE/(FCT(IA)*FCT(IB)*FCT(IC)*FCT(ID)*FCT(IE)*FCT(IF))
429 Z=Z+2
430 IF(Z.LE.ZMAX) GO TO 2
431 FASE=-DSIGN(1.0D0,CC)
432 IF(MOD(J1-J2,4).EQ.0) FASE=-FASE
433 CC=FASE*DSQRT(Y*FCT(JB)*FCT(JC)*FCT(JE)*CC*FCT(JA)*FCT(JD)*FCT(JE)
434 1*CC*FCT(JF)*FCT(JG)*FCT(JH)/FCT(JJ))
435 1 TRIXJ=CC
436 RETURN
437 3 PRINT 4
438 4 FORMAT('0 ERROR IN ARGUMENT OF TRIXJ')
439 CALL EXIT
440 5 PRINT 6
441
442
443 6 FORMAT('0 ERROR FACTORIALS EXCEEDED IN TIXJ')
444 CALL EXIT
445 END
446 C Seite 28
447
448
449 SUBROUTINE VPSI(RAD,THETA,FHI,NP,NB,M,ND,PSIR,PSITH,PSIFH,PN,U,V)
450 DIMENSION PSIR(1),PSITH(1),PSIFH(1),PN(1),U(1),V(1)
451 DOUBLE PRECISION RAD,R,THETA,T,FHI,F,PN,U,V,FCT,FR,FR1,FR2,B,C
452 COMPLEX*16 FC,FC1,FC2,PSIR,PSITH,PSIFH
453 COMMON/FAC/FCT(100),NRISK,NTRIX
454 R = RAD
455 T = THETA
456 F = FHI
457 NP1 = NP-1
458 K = M
459 N = ND
460 L = N/2
461 L1 = L+1
462 L2 = L+2
463 FC = (0.0D0,1.0D0)
464 IF(NB.EQ.1) GO TO 5
465 CALL BN(R,L1,U,V)
466 B = DABS(R*R*(U(2)*V(1)-U(1)*V(2))-1.0D0)
467 C = DABS(R*R*(U(L1)*V(L)-U(L)*V(L1))-1.0D0)
468 PRINT 3
469 3 FORMAT('0 BESSEL- NEUMAN- TEST FOR VPSI')
470 PRINT 4,B,C
471 4 FORMAT(2D25.15)
472 5 CONTINUE
473 CALL LEG(T,K,L2,PN)
474 DO 42 I = 1,L
475 I1 = I+1
476 I2 = I+2
477 IF(I-K)6,7,7
478 6 FR = 0.0D0
479 GO TO 10
480 7 IF(K)8,8,9
481 8 FR = DSQRT(DFLOAT(2*I+1)/(DFLOAT(I*I1)*16.0D0*DATAN(1.0D0)))
482 GO TO 10
483 9 FR = DSQRT(DFLOAT(2*I+1)*FCT(I1-K)/(DFLOAT(I*I1)*FCT(I1+K)*8.0D0*
484 1DATAN(1.0D0)))
485 10 CONTINUE
486 IF(T)16,16,19
487 16 CONTINUE
488 IF(K-1)18,17,18
489 17 CONTINUE
490 FR1 = DSQRT(DFLOAT(2*I+1)/(32.0D0*DATAN(1.0D0)))
491 FR2 = DSQRT(DFLOAT(2*I+1)/(32.0D0*DATAN(1.0D0)))
492 GO TO 20
493 18 CONTINUE
494 FR1 = 0.0D0
495 FR2 = 0.0D0
496 GO TO 20
497 19 CONTINUE
498 FR1 = FR*PN(I1)*DFLOAT(K)/DSIN(T)
499 FR2 = -FR*(DFLOAT(I1)*DCOS(T)*PN(I1)-DFLOAT(I1-K)*PN(I2))/DSIN(T)
500 20 CONTINUE
501 GO TO (30,31,32),NB
502 30 CONTINUE
503 FC1 = (-FC)**I1
504 FC2 = (-FC)**I
505 GO TO 33
506 31 CONTINUE
507 C Seite 29
508
509
510 FC1 = DCMPLX(U(I1),0.0D0)
511 FC2 = DCMPLX(U(I)-DFLOAT(I)*U(I1)/R,0.0D0)
512
513
514 GO TO 33
515 32 CONTINUE
516 FC1 = DCMPLX(U(I1),V(I1))
517 FC2 = DCMPLX(U(I)-DFLOAT(I)*U(I1)/R,V(I)-DFLOAT(I)*V(I1)/R)
518 33 CONTINUE
519 PSIR(2*I-1) = (0.0D0,0.0D0)
520 IF(NB-1)34,34,35
521 34 CONTINUE
522 PSIR(2*I) = (0.0D0,0.0D0)
523 35 CONTINUE
524 GO TO (36,39),NP1
525 36 CONTINUE
526 IF(NB-1)38,38,37
527 37 CONTINUE
528 PSIR(2*I) = FC1*DCMPLX(DFLOAT(I*I1)*FR*PN(I1)*DSIN(DFLOAT(K)*F)/R,
529 10.0D0)
530 38 CONTINUE
531 PSITH(2*I-1) = -FC1*DCMPLX(FR1*DSIN(DFLOAT(K)*F),0.0D0)
532 PSITH(2*I) = FC2*DCMPLX(FR2*DSIN(DFLOAT(K)*F),0.0D0)
533 PSIFH(2*I-1) = -FC1*DCMPLX(FR2*DCOS(DFLOAT(K)*F),0.0D0)
534 PSIFH(2*I) = FC2*DCMPLX(FR1*DCOS(DFLOAT(K)*F),0.0D0)
535 GO TO 42
536 39 CONTINUE
537 IF(NB-1)41,41,40
538 40 CONTINUE
539 PSIR(2*I) = FC1*DCMPLX(DFLOAT(I*I1)*FR*PN(I1)*DCOS(DFLOAT(K)*F)/R,
540 10.0D0)
541 41 CONTINUE
542 PSITH(2*I-1) = FC1*DCMPLX(FR1*DCOS(DFLOAT(K)*F),0.0D0)
543 PSITH(2*I) = FC2*DCMPLX(FR2*DCOS(DFLOAT(K)*F),0.0D0)
544 PSIFH(2*I-1) = -FC1*DCMPLX(FR2*DSIN(DFLOAT(K)*F),0.0D0)
545 PSIFH(2*I) = -FC2*DCMPLX(FR1*DSIN(DFLOAT(K)*F),0.0D0)
546 42 CONTINUE
547 RETURN
548 END
549 C Seite 30
550
551
552 SUBROUTINE VKOEF (BHETA,NP,M,ND,AP,PN)
553 DIMENSION PN(1),AP(1)
554 DOUBLE PRECISION BHETA,U,DI,PN,FR,FCT
555 COMPLEX*16 FC1,FC2,fc3,AP
556 COMMON/FAC/FCT(100),NRISK,NTRIX
557 N = NP
558 N1 = N+1
559 K = M
560 L = ND/2
561 L2 = L+2
562 U = BHETA
563 DI = DATAN(1.0D0)
564 FC1 = (0.0D0,1.0D0)
565 CALL LEG(U,K,L2,PN)
566 IF(U)10,10,20
567 10 DO 100 I = 1,L
568 IF(K-1)11,12,11
569 11 AP(2*I-1) = (0.0D0,0.0D0)
570 AP(2*I) = (0.0D0,0.0D0)
571 GO TO 100
572 12 FR = DSQRT(8.0D0*DI*DFLOAT(2*I+1))
573 FC2 = DCMPLX(FR,0.0D0)
574 GO TO (13,14),N1
575 13 AP(2*I-1) = -(FC1**I)*FC2
576 AP(2*I) = -(FC1**(I+1))*FC2
577 GO TO 100
578 14 AP(2*I-1) = (FC1**I)*FC2
579 AP(2*I) = (FC1**(I+3))*FC2
580 100 CONTINUE
581 GO TO 300
582 20 DO 200 I = 1,L
583 I1 = I+1
584 I2 = I+2
585 IF(K)21,21,30
586 21 FR = 4.0D0*DSQRT((DI*DFLOAT((2*I+1)*(I+1)))/DFLOAT(I))*(DCOS(U)*
587 1PN(I1)-PN(I2))/DSIN(U)
588 FC2 = DCMPLX(FR,0.0D0)
589 GO TO (22,23),N1
590 22 AP(2*I-1) = (FC1**I)*FC2
591 AP(2*I) = (0.0D0,0.0D0)
592 GO TO 200
593 23 AP(2*I-1) = (0.0D0,0.0D0)
594 AP(2*I) = (FC1**(I+1))*FC2
595 GO TO 200
596 30 IF(I-K)34,31,31
597 31 FR = 4.0D0*DSQRT((2.0D0*DI*DFLOAT(2*I+1)*FCT(I-K+1))/(DFLOAT(I*(I
598 1+1))*FCT(I+K+1)))
599 FC2 = DCMPLX((FR*(DFLOAT(I+1)*DCOS(U)*PN(I1)-DFLOAT(I-K+1)*PN(I2)
600 1))/DSIN(U),0.0D0)
601 FC3 = DCMPLX((DFLOAT(K)*FR*PN(I1))/DSIN(U),0.0D0)
602 GO TO (32,33),N1
603 32 AP(2*I-1) = (FC1**I)*FC2
604 AP(2*I) = -(FC1**(I+1))*FC3
605 GO TO 200
606 33 AP(2*I-1) = (FC1**I)*FC3
607 AP(2*I) = (FC1**(I+1))*FC2
608 GO TO 200
609 34 AP(2*I-1) = (0.0D0,0.0D0)
610 AP(2*I) = (0.0D0,0.0D0)
611 200 CONTINUE
612
613
614 300 RETURN
615 END
616 C Seite 31
617
618
619 SUBROUTINE VR(AR,NP,M,ND,R1,R2,R3,R4,X1,X2,Y1,Y2)
620 DIMENSION R1(ND,1),R2(ND,1),R3(ND,1),R4(ND,1),X1(1),X2(1),Y1(1),Y2
621 1(1)
622 DOUBLE PRECISION AR,U,V,F3J1,F3J2,F3J3,FACT1,FACT2,FACT3,FACT4,
623 1FACT5,FACT6,FACT7,FACT8,FACT9,FACT10,ACR1,ACR2,ACR3,ACR4,X1,X2,
624 1Y1,Y2,B1,B2,CN1,CN2,FCT,trixj
625 COMPLEX*16 R1,R2,R3,R4,W1,W2,W3,W4,W5,W6,W7,W8
626 COMMON/FAC/FCT(100),NRISK,NTRIX
627 NP1 = NP+1
628 N = M
629 NK = ND+1
630 U = AR
631 V = 0.5D0*AR
632 CALL BN(U,NK,X1,Y1)
633 CALL BN(V,NK,X2,Y2)
634 L1 = NK
635 L = NK-1
636 B1 = DABS(U**2*(X1(2)*Y1(1)-X1(1)*Y1(2))-1.0D0)
637 B2 = DABS(V**2*(X1(2)*Y2(1)-X2(1)*Y2(2))-1.0D0)
638 CN1 = DABS(U**2*(X1(L1)*Y1(L)-X1(L)*Y1(L1))-1.0D0)
639 CN2 = DABS(V**2*(X2(L1)*Y2(L)-X2(L)*Y2(L1))-1.0D0)
640 PRINT 89
641 89 FORMAT('0 BESSEL- NEUMAN- TEST FOR VR')
642 PRINT 90, B1,CN1,B2,CN2
643 90 FORMAT(4D25.15)
644 DO 7 I = 1,L
645 DO 7 J = 1,L
646 R1(I,J) = (0.0D0,0.0D0)
647 R2(I,J) = (0.0D0,0.0D0)
648 R3(I,J) = (0.0D0,0.0D0)
649 R4(I,J) = (0.0D0,0.0D0)
650 7 CONTINUE
651 NR = L/2
652 DO 200 I = 1,NR
653 DO 200 J = 1,NR
654 IF(N-I)8,8,200
655 8 IF(N-J)9,9,200
656 9 L1 = I+J+1
657 W1 = (0.0D0,0.0D0)
658 W2 = (0.0D0,0.0D0)
659 W3 = (0.0D0,0.0D0)
660 W4 = (0.0D0,0.0D0)
661 W5 = (0.0D0,0.0D0)
662 W6 = (0.0D0,0.0D0)
663 W7 = (0.0D0,0.0D0)
664 W8 = (0.0D0,0.0D0)
665 DO 100 L = 1,L1
666 K = L-1
667 IF(K-IABS(I-J))100,10,10
668 10 CONTINUE
669 IF(N)11,11,12
670 11 IF(MOD((I+J+K),2).NE.0) GO TO 100
671 FACT1 = 1.0D0
672 GO TO 13
673 12 FACT1 = DFLOAT((-1)**N)
674 13 J1 = 2*I
675 J2 = 2*J
676 J3 = 2*K
677 M1 = 2*N
678 C Seite 32
679
680
681 F3J1 = TRIXJ(J1,J2,J3,M1,FCT,NRISK,NTRIX)
682 FACT2 = 0.5D0*DFLOAT(2*K+1)*F3J1*
683
684
685 1DSQRT(DFLOAT((2*I+1)*(2*J+1))/DFLOAT(I*(I+1)*J*(J+1)))
686 FACT3 = FACT1*FACT2
687 IF(MOD((J-I+K),2).NE.0) GO TO 27
688 IF(J-I+K)25,23,24
689 23 FACT4 = 1.0D0
690 GO TO 26
691 24 FACT4 = DFLOAT((-1)**((J-I+K)/2))
692 GO TO 26
693 25 FACT4 = DFLOAT((-1)**((I-J-K)/2))
694 26 J1 = 2*I
695 J2 = 2*J
696 J3 = 2*K
697 M1 = 0
698 F3J2 = TRIXJ(J1,J2,J3,M1,FCT,NRISK,NTRIX)
699 FACT5 = DFLOAT(I*(I+1)+J*(J+1)-K*(K+1))*F3J2
700 GO TO 28
701 27 FACT6 = 0.0D0
702 GO TO 29
703 28 FACT6 = FACT4*FACT5
704 29 IF(K-IABS(I-J)-1)44,30,30
705 30 IF(MOD((J-I+K+1),2).NE.0) GO TO 44
706 IF(J-I+K+1)33,31,32
707 31 FACT7 = 1.0D0
708 GO TO 41
709 32 FACT7 = DFLOAT((-1)**((J-I+K+1)/2))
710 GO TO 41
711 33 FACT7 = DFLOAT((-1)**((I-J-K-1)/2))
712 41 J1 = 2*I
713 J2 = 2*J
714 J3 = 2*(K-1)
715 M1 = 0
716 F3J3 = TRIXJ(J1,J2,J3,M1,FCT,NRSIK,NTRIX)
717 IF(IABS(I-J))42,42,43
718 42 FACT8 = DFLOAT(K)*DSQRT(DFLOAT((I+J-1)**2-K**2))*F3J3
719 GO TO 45
720 43 FACT8 = DSQRT(DFLOAT((K**2-IABS(I-J)**2)*((I+J+1)**2-K**2)))*F3J3
721 GO TO 45
722 44 FACT9 = 0.0D0
723 GO TO 50
724 45 FACT9 = FACT7*FACT8
725 50 CONTINUE
726 IF(K)51,51,52
727 51 FACT10 = 1.0D0
728 GO TO 53
729 52 FACT10 = DFLOAT((-1)**K)
730 53 ACR1 = FACT3*FACT6
731 ACR2 = FACT10*ACR1
732 ACR3 = FACT3*FACT9
733 ACR4 = FACT10*ACR3
734 K1 = K+1
735 W1 = W1+DCMPLX(ACR1*X1(K1),0.0D0)
736 W2 = W2+DCMPLX(ACR1*X1(K1),ACR1*Y1(K1))
737 W3 = W3+DCMPLX(ACR2*X1(K1),ACR2*Y1(K1))
738 W4 = W4+DCMPLX(ACR1*X2(K1),0.0D0)
739 IF(N)100,100,99
740 C Seite 33
741
742
743 99 CONTINUE
744 W5 = W5+DCMPLX(ACR3*X1(K1),0.0D0)
745 W6 = W6+DCMPLX(ACR3*X1(K1),ACR3*Y1(K1))
746 W7 = W7+DCMPLX(ACR4*X1(K1),ACR4*Y1(K1))
747
748
749 W8 = W8+DCMPLX(ACR3*X2(K1),0.0D0)
750 100 CONTINUE
751 R1(2*I-1,2*J-1) = W1
752 R2(2*I-1,2*J-1) = W2
753 R3(2*I-1,2*J-1) = W3
754 R4(2*I-1,2*J-1) = W4
755 R1(2*I,2*J) = W1
756 R2(2*I,2*J) = W2
757 R3(2*I,2*J) = W3
758 R4(2*I,2*J) = W4
759 GO TO (101,103),NP1
760 101 CONTINUE
761 IF(N)200,200,102
762 102 CONTINUE
763 R1(2*I-1,2*J) = W5
764 R2(2*I-1,2*J) = W6
765 R3(2*I-1,2*J) = W7
766 R4(2*I-1,2*J) = W8
767 R1(2*I,2*J-1) = -W5
768 R2(2*I,2*J-1) = -W6
769 R3(2*I,2*J-1) = -W7
770 R4(2*I,2*J-1) = -W8
771 GO TO 200
772 103 CONTINUE
773 IF(N)200,200,104
774 104 CONTINUE
775 R1(2*I-1,2*J) = -W5
776 R2(2*I-1,2*J) = -W6
777 R3(2*I-1,2*J) = -W7
778 R4(2*I-1,2*J) = -W8
779 R1(2*I,2*J-1) = W5
780 R2(2*I,2*J-1) = W6
781 R3(2*I,2*J-1) = W7
782 R4(2*I,2*J-1) = W5
783 200 CONTINUE
784 RETURN
785 END
786 C Seite 34
787
788
789 SUBROUTINE MCNV(A,N,D,L,M)
790 DIMENSION A(1),L(1),M(1)
791 COMPLEX*16 A,D,BIGA,HOLD
792 D = (1.0D0,0.0D0)
793 NK=-N
794 DO 80 K=1,N
795 NK=NK+N
796 L(K)=K
797 M(K)=K
798 KK=NK+K
799 BIGA=A(KK)
800 DO 20 J=K,N
801 IZ=N*(J-1)
802 DO 20 I=K,N
803 IJ=IZ+1
804 10 IF(CDABS(BIGA)-CDABS(A(IJ))) 15,20,20
805 15 BIGA=A(IJ)
806 L(K)=I
807 M(K)=J
808 20 CONTINUE
809 J=L(K)
810 IF(J-K) 35,35,25
811 25 KI=K-N
812 DO 30 I=1,N
813 KI=KI+N
814 HOLD=-A(KI)
815 JI=KI-K+J
816 A(KI)=A(JI)
817 30 A(JI) =HOLD
818 35 I=M(K)
819 IF(I-K) 45,45,38
820 38 JP=N*(I-1)
821 DO 40 J=1,N
822 JK=NK+J
823 JI=JP+J
824 HOLD=-A(JK)
825 A(JK)=A(JI)
826 40 A(JI) =HOLD
827 45 IF(CDABS(BIGA)) 48,46,48
828 46 D = (0.0D0,0.0D0)
829 RETURN
830 48 DO 55 I=1,N
831 IF(I-K) 50,55,50
832 50 IK=NK+1
833 A(IK)=A(IK)/(-BIGA)
834 55 CONTINUE
835 DO 65 I=1,N
836 IK=NK+1
837 HOLD=A(IK)
838 IJ=I-N
839 DO 65 J=1,N
840 IJ=IJ+N
841 IF(I-K) 60,65,60
842 60 IF(J-K) 62,65,62
843 62 KJ=IJ-I+K
844 A(IJ)=HOLD*A(KJ)+A(IJ)
845 65 CONTINUE
846 KJ=K-N
847 DO 75 J=1,N
848 C Seite 35
849
850
851 KJ=KJ+N
852
853
854 IF(J-K) 70,75,70
855 70 A(KJ)=A(KJ)/BIGA
856 75 CONTINUE
857 D=D*BIGA
858 A(KK) = (1.0D0,0.0D0)/BIGA
859 80 CONTINUE
860 K=N
861 100 K=(K-1)
862 IF(K) 150,150,105
863 105 I=L(K)
864
865
866 IF(I-K) 120,120,108
867 108 JQ=N*(K-1)
868 JR=N*(I-1)
869 DO 110 J=1,N
870 JK=JQ+J
871 HOLD=A(JK)
872 JI=JR+J
873 A(JK)=-A(JI)
874 110 A(JI) =HOLD
875 120 J=M(K)
876 IF(J-K) 100,100,125
877 125 KI=K-N
878 DO 130 I=1,N
879 KI=KI+N
880 HOLD=A(KI)
881 JI=KI-K+J
882 A(KI)=-A(JI)
883 130 A(JI) =HOLD
884 GO TO 100
885 150 RETURN
886 END
887 C Seite 36
888
889
890 SUBROUTINE COND(M,ND,RE,IM)
891 DIMENSION RE(ND,1),IM(ND,1)
892 DOUBLE PRECISION RE,IM,SCALE
893 MD = 1
894 IF(M.GT.1) MD = 2*M-1
895 MD1 = MD+1
896 NBGR = ND
897 NROW = NBGR
898 DO 60 KR = MD1,NBGR
899 SCALE = 1.0D0/IM(NROW,NROW)
900 DO 8 LC = MD,NBGR
901 RE(NROW,LC) = SCALE*RE(NROW,LC)
902 IM(NROW,LC) = SCALE*IM(NROW,LC)
903 8 CONTINUE
904 MROW = NROW-1
905 DO 20 MR = MD,MROW
906 SCALE = IM(MR,NROW)
907 DO 16 MC = MD,NBGR
908 RE(MR,MC) = RE(MR,MC)-SCALE*RE(NROW,MC)
909 IM(MR,MC) = IM(MR,MC)-SCALE*IM(NROW,MC)
910 16 CONTINUE
911 20 CONTINUE
912 NROW = NROW-1
913 60 CONTINUE
914 NROW = NBGR-1
915 DO 80 I = 1,NROW
916 IB = I+1
917 DO 72 J = IB,NBGR
918 IM(I,J) = 0.0D0
919 72 CONTINUE
920 80 CONTINUE
921 RETURN
922 END
923 C Seite 37
924
925
926 SUBROUTINE ORTHO(M,ND,RE,IM,X,Y)
927 DIMENSION RE(ND,1),IM(ND,1),X(1),Y(1)
928 DOUBLE PRECISION RE,IM,X,Y,SUM1,SUM2
929 MD = 1
930 IF(M.GT.1) MD = 2*M-1
931 NBGR = ND
932 SUM1 = 0.0D0
933 DO 20 K = MD,NBGR
934 SUM1 = SUM1+RE(NBGR,K)**2+IM(NBGR,K)**2
935 20 CONTINUE
936 SUM1 = DSQRT(SUM1)
937 DO 28 K = MD,NBGR
938 RE(NBGR,K) = RE(NBGR,K)/SUM1
939 IM(NBGR,K) = IM(NBGR,K)/SUM1
940 28 CONTINUE
941 NMI = NBGR-1
942 NROW = NBGR
943 DO 100 I = MD,NMI
944 NROW = NROW-1
945 MROW = NROW
946 DO 36 K = MD,NBGR
947 X(K) = RE(NROW,K)
948 Y(K) = IM(NROW,K)
949 36 CONTINUE
950 DO 80 J = NROW,NMI
951 SUM1 = 0.0D0
952 SUM2 = 0.0D0
953 MROW = MROW+1
954 DO 40 K = MD,NBGR
955 SUM1 = SUM1+RE(MROW,K)*RE(NROW,K)+IM(MROW,K)*IM(NROW,K)
956 SUM2 = SUM2+RE(MROW,K)*IM(NROW,K)-IM(MROW,K)*RE(NROW,K)
957 40 CONTINUE
958 DO 48 K = MD,NBGR
959 X(K) = X(K)-SUM1*RE(MROW,K)+SUM2*IM(MROW,K)
960 Y(K) = Y(K)-SUM1*IM(MROW,K)-SUM2*RE(MROW,K)
961 48 CONTINUE
962 80 CONTINUE
963 SUM1 = 0.0D0
964 DO 84 K = MD,NBGR
965 SUM1 = SUM1+X(K)**2+Y(K)**2
966 84 CONTINUE
967 SUM1 = DSQRT(SUM1)
968 DO 88 K = MD,NBGR
969 RE(NROW,K) = X(K)/SUM1
970 IM(NROW,K) = Y(K)/SUM1
971 88 CONTINUE
972 100 CONTINUE
973 RETURN
974 END
975 C Seite 38
976
977
978 SUBROUTINE PERFT(NP,M,NR,ND,AR,AI,BR,BI,T,RE,IM,X,Y)
979 DIMENSION AR(NR,1),AI(NR,1),BR(NR,1),BI(NR,1),T(ND,1),RE(ND,1),
980 1IM(ND,1),X(1),Y(1)
981 DOUBLE PRECISION AR,AI,BR,BI,RE,IM,X,Y,FAC
982 COMPLEX*16 T
983 MD = 1
984 IF(M.GT.1) MD = 2*M-1
985 NR = ND/2
986 FAC = 1.0D0
987 IF(NP.GT.0) FAC = -1.0D0
988 DO 10 I = 1,NR
989 DO 10 J = 1,NR
990 RE(2*I-1,2*J-1) = AR(I,J)
991 RE(2*I-1,2*J) = FAC*BR(I,J)
992 RE(2*I,2*J-1) = FAC*BR(I,J)
993 RE(2*I,2*J) = -AR(I,J)
994 IM(2*I-1,2*J-1) = AI(I,J)
995 IM(2*I-1,2*J) = FAC*BI(I,J)
996 IM(2*I,2*J-1) = FAC*BI(I,J)
997 IM(2*I,2*J) = -AI(I,J)
998 IF(1.EQ.J) IM(2*I,2*I) = 1.0D0-AI(I,I)
999 10 CONTINUE
1000 CALL COND(M,ND,RE,IM)
1001 CALL ORTHO(M,ND,RE,IM,X,Y)
1002 DO 11 I = 1,ND
1003 DO 11 J = 1,ND
1004 T(I,J) = (0.0D0,0.0D0)
1005 11 CONTINUE
1006 DO 12 I = MD,ND
1007 DO 12 J = MD,ND
1008 DO 12 K = MD,ND
1009 T(I,J) = T(I,J)-DCMPLX(RE(K,J)*RE(K,J),-IM(K,I)*RE(K,J))
1010 12 CONTINUE
1011 RETURN
1012 END
1013 C Seite 39
1014
1015
1016 SUBROUTINE LINE(THETA,NIP,C,ALPHA,R,DR)
1017 DOUBLE PRECISION THETA,C,ALPHA,R,DR
1018 R = C*DSIN(ALPHA)/DSIN(THETA+DFLOAT(NIP)*ALPHA)
1019 DR = -R/DTAN(THETA+DFLOAT(NIP)*ALPHA)
1020 RETURN
1021 END
1022
1023
1024 SUBROUTINE TRCIR(STH,CTH,C,A,R,DR)
1025 DOUBLE PRECISION STH,CTH,C,A,R,DR,X,ST,CT
1026 ST = C*STH
1027 CT = C*CTH
1028 X = DSQRT(A**2-ST**2)
1029 R = CT+X
1030 DR = -ST-ST*CT/X
1031 RETURN
1032 END
1033
1034
1035 SUBROUTINE ELLIPS(STH,CTH,AZ,BRA,R,DR)
1036 DOUBLE PRECISION STH,CTH,AZ,BRA,R,DR
1037 R = AZ/DSQRT (CTH**2+(AZ*STH/BRA)**2)
1038 DR = R**3*STH*CTH*(1.0D0-(AZ/BRA)**2)/AZ**2
1039 RETURN
1040 END
1041
1042
1043 SUBROUTINE TRELLI(STH,CTH,AZ,BRA,C,R,DR)
1044 DOUBLE PRECISION STH,CTH,AZ,BRA,C,R,DR,NE,RO
1045 NE = (AZ*STH)**2+(BRA*CTH)**2
1046 RO = NE-(C*STH)**2
1047 RO = DSQRT(RO)
1048 R = (BRA**2*C*CTH+AZ*BRA*RO)/NE
1049 DR = -2.0D0*STH*CTH*(AZ**2-BRA**2)*R/NE+(-BRA**2*C*STH+STH*CTH*AZ*
1050 1BRA*(AZ**2-BRA**2-C**2)/RO)/NE
1051 RETURN
1052 END